/* 
*  dpmta_legendre.h - inline legendre calcuation macros
*
*  w.t.rankin
*
*  macros taken from:
*
*     PMTAlegendre.h
*        Version 3.0, February 20, 1994
*
*  Copyright (c) 1994 Duke University
*  All rights researved
*
*  RCS Info: $Id: dpmta_legendre.h,v 1.1.1.1 2004/12/17 13:19:16 jspturk Exp $
*
*/

/*
 * revision history:
 *
 * $Log: dpmta_legendre.h,v $
 * Revision 1.1.1.1  2004/12/17 13:19:16  jspturk
 * adding particle in LEAVES with link list instead of an array.
 *
 * Revision 1.2  1994/10/14  04:42:14  wrankin
 * modifications to support new multipole library
 *
 * Revision 1.1  1994/10/09  18:25:37  wrankin
 * Initial revision
 *
 *
 */



/****************************************************************
*
*  Legendre() sets up the Legendre Polynomial
*  adapted from Numerical Recipes in C
* 
*  P is type Mtype, xval is type double, and p is type unsigned int
*  p is the number of terms in the multipole expansion
*
*/

#define Legendre(P,p,xval)                                              \
{                                                                       \
    int Li, Lj;                                                         \
    double negterm, oddfact, nextodd, sqroot, sqrtterm;                 \
                                                                        \
    negterm = 1.0;                                                      \
    oddfact = 1.0;                                                      \
    nextodd = 1.0;                                                      \
    sqroot = sqrt(1.0 - xval*xval);                                     \
    sqrtterm = 1.0;                                                     \
    for(Li=0;Li < p;Li++){                                              \
        P[Li][Li] = negterm*oddfact*sqrtterm;                           \
        negterm *= -1.0;                                                \
        oddfact *= nextodd;                                             \
        nextodd += 2.0;                                                 \
        sqrtterm *= sqroot;                                             \
        if(Li < p-1){                                                     \
        P[Li+1][Li] = xval * (double)(2*Li+1) * P[Li][Li];              \
            for(Lj=Li+2;Lj < p;Lj++){                                  \
                P[Lj][Li] = (xval*(double)(2*Lj-1)*P[Lj-1][Li] -        \
                    (double)(Lj+Li-1)*P[Lj-2][Li])/(double)(Lj-Li);     \
            }                                                           \
        }                                                               \
    }                                                                   \
}

